package fem2.element;

import inf.jlinalg.ArrayVector;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IMatrixReadOnly;
import inf.jlinalg.IVector;
import inf.jlinalg.IVectorReadOnly;
import inf.math.LinAlg;
import fem2.Edge;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.enu.IntegrationType;
import fem2.enu.MassMatrixType;

/**
 * TODO: element is under construction
 * 
 * @author hbui
 * 
 */
public class NonlinearTrusssElement extends StructuralElement {

	/**
	 * @param mp
	 *            mesh part
	 * @param mm
	 *            material model
	 * @param preStress
	 *            pre-stress in the truss
	 */
	public NonlinearTrusssElement(MeshPart mp, MaterialModel mm, double preStress) {
		super(mp, mm);
		if (!(mp instanceof Edge)) {
			throw new Error("truss element must be defined on an edge");
		}
	}

	/**
	 * @return the mesh part as Edge
	 */
	public Edge getEdge() {
		return (Edge) getMeshPart();
	}

	@Override
	public void computeB(MeshPart mp, IMatrix b, double... xi) {
		int n = mp.getBasis().getDimension();
		Edge e = getEdge();
		double[][] gGradients = e.gGradientsAt(xi);
		double[][] ju = mp.gJacUAt(xi);
		double u11 = ju[0][0];
		for (int i = 0; i < n; i++) {
			double g1 = gGradients[i][0];
			b.set(0, i, g1 + g1 * u11);
		}
	}

	@Override
	public IMatrixReadOnly computeKte(IMatrix ke) {
		// TODO Auto-generated method stub
		return super.computeKte(ke);
	}

	@Override
	public IVectorReadOnly computeRie(IVector rie) {
		// TODO Auto-generated method stub
		return super.computeRie(rie);
	}

	@Override
	public IVectorReadOnly computeRes(IVector res) {
		// TODO Auto-generated method stub
		return super.computeRes(res);
	}

	@Override
	public IMatrixReadOnly computeMe(IMatrix me, MassMatrixType massType) {
		// TODO Auto-generated method stub
		return super.computeMe(me, massType);
	}

	@Override
	public IMatrixReadOnly computeKge(IMatrix ke) {
		MeshPart mp = getMeshPart();
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		int ng = getNumIntegrationPoints();
		double a = getMaterialModel().getSectionProperty();
		double s = 1.0;

		if (getIntegrationType() == IntegrationType.SELECTIVE_REDUCED_INTEGRATION) {
			s = 0.5; // to account for more Gauss points in the bulk
		}

		for (int i = 0; i < ng; i++) {
			double[] xi = getIntegrationPoint(i);
			double w = getIntegrationWeight(i);

			double[][] jac = mp.jacXAt(xi);
			double detJ = LinAlg.norm2(dim, jac[0]);

			double[][] gv = mp.gGradientsAt(xi);
			IVector sigma = getStress(i);

			for (int j = 0; j < n; j++) {
				for (int k = 0; k < n; k++) {

					double ni1 = gv[j][0];
					double nj1 = gv[k][0];

					double f = sigma.get(0) * ni1 * nj1;
					double coeff = w * detJ * f * a * s;
					for (int l = 0; l < dim; l++) {
						ke.add(dim * j + l, dim * k + l, coeff);
					}
				}
			}
		}
		return ke;
	}

	@Override
	public IVector computeRegularStrain(double... xi) {
		MeshPart mp = getMeshPart();
		double[][] ju = mp.gJacUAt(xi);
		int ntrs = 1;
		IVector e = new ArrayVector(ntrs);
		double u11 = ju[0][0];
		e.set(0, 0.5 * (2 * u11 + u11 * u11));
		return e;
	}

	@Override
	public IVector computeStrain(double... xi) {
		return computeRegularStrain(xi);
	}

	// @Override
	// public double computeStrainEnergy() {
	// /*
	// * TODO
	// */
	// throw new NotImplementedError();
	// }

}
